;+
; NAME: 
;       H2_DENSITY
;
; PURPOSE: 
;       Return the density of H2 in the Milky Way as a function of R,Z
;       in cylindrical coordinates.  Uses a new double-gaussian model
;       fit to the Bronfman (1988) CO data for the Northern Galactic
;       Plane.  NOTE: All values scaled to R0 of Reid (2009).
;
; CALLING SEQUENCE:  
;       rho = H2_DENSITY(r,z)
;
; INPUTS:
;       R, Z -- Galactocentric radius and height above/below the midplane
;               in pc.
;
; KEYWORD PARAMETERS:
;       NONE
;
; OUTPUTS:
;       RHO -- volume mass density of H2
;
; MODIFICATION HISTORY:
;
;      Modified: 09/07/10, TPEB -- Renamed the routine, and made all
;                                  galactocentric distances scalable
;                                  by the current version of R0.
;      Modified: 09/27/10, TPEB -- Changed the model fit to a
;                                  double-gaussian fit to just the
;                                  Northern plane of Bronfman 1988.
;
;-

FUNCTION H2_DENSITY, R, Z
  
  COMPILE_OPT OBSOLETE
  
  ;; Get galactic params
  defsysv, '!MW', exists = exists
  IF NOT exists THEN galactic_params 
  
  R0 = !mw.r0
  
  
  
  
  ;; Fiddle-able parameters
  ;; Bronfmann et al (2000) Gaussian Parameters (Table 5) -- North
  maxden        = 5.5d                      ;; M_sun / pc^2   \pm 0.4
  den_sc        = 1.4d                      ;; M_sun / pc^2   at Solar Circle
  g_cen         = 0.57d * R0                ;; \pm 0.01
  g_fwhm        = 0.37d * R0                ;; \pm 0.03 
  out_sl        = 0.29d * R0                ;; \pm 0.02
  inner_g_hwhm  = 59.d                      ;; pc
  outer_g_scale = 6.7d3 * R0 / 8.5d3        ;; Scaled value to new R0
  trans         = 0.695d * R0 ;; pc  Transition from inner model to outer 
  
  ;; Scale Height
  scaleht = 2*inner_g_hwhm/sqrt(8*alog(2)) * (1 > exp((r-R0)/(outer_g_scale)) )
  
  ;; Case 1 -- Within Transition Point
  sd1 = maxden * exp(- (r-g_cen)^2 / (2*g_fwhm^2/(8*alog(2))))
  rho1 = sd1/sqrt(2*!pi)/scaleht*exp(-z^2/(2*scaleht^2))
  
  ;; Case 2 -- Outside Transition Point
  sd2 = den_sc * exp(-(r-R0)/out_sl)
  rho2 =  sd2/sqrt(2*!pi)/scaleht*exp(-z^2/(2*scaleht^2))
  
  rho = (r le trans)*rho1+(r gt trans)*rho2
  
  RETURN, rho
END
